function error_value = process_error_2d(node, elem, basis_type, pde, U, error_type)

gauss_weights_tri = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];
gauss_weights_rec = [0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.79012345679012345679012345679012, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864];

psi = reference_basis_points_2d(basis_type);
phi_x = reference_basis_points_2d(basis_type, "dx");
phi_y = reference_basis_points_2d(basis_type, "dy");
ref_trans_tri = reference_basis_points_2d("P1");
ref_trans_rec = reference_basis_points_2d("Q1");

switch error_type
    case "L2"
        sum_n = 0;
        switch basis_type
            case {"P1", "P1b", "P2"}
                for n = 1:size(elem,1)
                    E = node(elem(n,:),:);
                    V = E(1:3,1:2);
                    uh_E = zeros(1,9);
                    for beta = 1:size(elem,2)
                        uh_E = uh_E + U(elem(n,beta)).*psi{beta};
                    end
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    ref_points = [ref_trans_tri{1}(:), ref_trans_tri{2}(:), ref_trans_tri{3}(:)]*V;
                    f = (pde.u(ref_points(:,1),ref_points(:,2))' - uh_E).^2;
                    r = (gauss_weights_tri * f') .* abs(J_det);
                    sum_n = sum_n + r;
                end
            case {"Q1", "Q1b", "Q2"}
                V = E(1:4,1:2);
        end
        error_value = sqrt(sum_n);

    case "H1"
        sum_nx = 0;
        sum_ny = 0;
        switch basis_type
            case {"P1", "P1b", "P2"}
                for n = 1:size(elem,1)
                    E = node(elem(n,:),:);
                    V = E(1:3,1:2);
                    Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
                    J_det = det(Jacobi);
                    J_inv = inv(Jacobi);
                    uh_E = zeros(2,9);
                    phi = cellfun(@(x,y) J_inv' * [x; y], phi_x, phi_y, "UniformOutput", false);
                    for beta = 1:size(elem,2)
                        uh_E = uh_E + U(elem(n,beta)) .* phi{beta};
                    end
                    uh_Ex = uh_E(1,:);
                    uh_Ey = uh_E(2,:);
                    ref_points = [ref_trans_tri{1}(:), ref_trans_tri{2}(:), ref_trans_tri{3}(:)]*V;
                    f1 = (pde.dudx(ref_points(:,1),ref_points(:,2))' - uh_Ex).^2;
                    f2 = (pde.dudy(ref_points(:,1),ref_points(:,2))' - uh_Ey).^2;
                    rx = (gauss_weights_tri * f1') .* abs(J_det);
                    ry = (gauss_weights_tri * f2') .* abs(J_det);
                    sum_nx = sum_nx + rx;
                    sum_ny = sum_ny + ry;
                end
            case {"Q1", "Q1b", "Q2"}
                V = E(1:4,1:2);
        end
        error_value = sqrt(sum_nx + sum_ny);

    otherwise
        error("Invalid error type.");
end

end